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I. INTRODUCTION 



The second, third, and fourth moments of fluctuation velocity, 
temperature, and passive tracer concentrations, and their cross 
correlations have generally been used to estimate the magnitude of 
dispersive atmospheric transport. This applies to the transport of 
momentum, energy, humidity, and contaminants such as pollutants. 
To avoid overestimates of transport, we must also be able to 
distinguish between wave induced and turbulent contributions to 
these fluctuation quantities, since periodic wave motion is much 
less dispersive than aperiodic turbulence. This applies 
particularly in transitionally turbulent regimes and regions such 
as the nocturnal atmospheric boundary layer and convective 
entrainment zone, where the mean flow is punctuated by both 
intermittent turbulence and waves. However, current algorithms for 
separating waves from turbulence are either unreliable, not 
operationally useful, or both. a facile wave/turbulence 
discriminant has not been demonstrated in the literature. 

Here we study D^, a self-affine, multi-fractal which may be the 
first operationally useful wave/turbulence discriminant for typical 
time series data. We discuss and demonstrate its advantages over 
a self-similar fractal, Fourier transforms, and other standard 
turbulence measures, such as Richardson number (Rj) , Brunt-Vaisala 
frequency (BVF) , buoyancy length scale (1 b) » variances, and 
turbulent kinetic energy (TKE) . We also discuss this new method's 
operational advantages over the phase averaging technique. For the 
comparisons we use nocturnal data from the 300 meter tower at the 
Boulder Atmospheric Observatory. We also suggest that be tested 
as a general measure of the degree of chaos in systems ranging from 
waves and limit cycles to standard Lorenz, Henon, and Poincare 
maps, to transitional turbulence, and turbulence displaying an 
extended inertial subrange. 

Fractal self-similarity exists when a geometric feature re-appears 
at successively smaller scales, where each successive scale is 
separated by a constant scaling factor (Mandelbrot, 1977). In 
fluids where turbulence lacks rigid geometric structure, we must 
instead test the statistical properties describing the ensemble 
mean geometry of the flow structure for possible self-similarity. 

Previous studies have suggested strong potential for a fractal 
self-similarity approach to wave/turbulence discrimination. 
Feigenbaum (1980) discovered that low order turbulence, i.e., 
chaos, in viscousless systems possesses an underlying self- 
similarity and hence fractally intermittent character. He showed 
that this self-simi larity in the density of Lagrangian particle 
streamlines appears only when the attractor, the point or area 
about which the particle motion occurs in real or phase space 
becomes unstable, i.e., strange. For waves or even limit cycle 
motion about fixed attractors, self-similarity does not occur. 
Pertinent to time series analysis, this self-similarity appears in 
real as well as phase space. 
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with regard to real viscous fluids, Kamada (1988a, b) showed a 
length scale self-similarity in atmospheric transitional turbulence 
for the capping inversion above a convective mixed layer. 
Sreenivasan and Meneveau (1986), and Presad and Sreenivasan (1990) 
found that the interfacial convolutions between turbulent and non- 
turbulent regions of a shear flow are self-similar and thus have a 
fractal dimension. However, pure self-similarity does not apply 
across all atmospheric turbulence scales because eddies are 
increasingly compressed vertically with scale, due to stability. 
However, Lovejoy and Schertzer (1984) included this compression by 
using a multi-fractal, self-affine stretch modification to the 
standard self-similar fractal. Packard et al. (1980) and Pawelzik 
and Schuster (1987) also studied the fractal dimensions of 
attractors inferred from a time series. 

However, the aim here is to find and utilize the fractal dimension 
of the time series itself, not the dimension of the phase space 
attractor or the dimension of the turbulence field in real space. 
Since a time series is a digital sampling, a time series trace of 
a fractal process should also show fractal characteristics. For 
example. Carter et al. (1986) found that the known fractal 
characteristics of cloud geometry in real space were also evident 
in the time series of their infrared intensity versus azimuth 
viewing angle. 

In real applications fractal analysis has potential advantages over 
standard spectral methods. 1) Less data manipulation is required. 
Standard Fourier transform techniques require that the data be 
periodic within the data window to avoid the introduction of 
fictitious high frequency noise. This requires tapering the data 
within the window so that its endpoints have the same value. This 
is not required for fractal analysis. 2) Fractal analysis' 
biggest advantage is that sudden amplitude changes, involving such 
phenomena as truncated gravity wave trains, breaking Kelvin- 
Helmholtz instabilities, or turbulent ramp structures pose no 
problems. In Fourier analyses a linear combination of wave numbers 
much higher than the time resolution is required to accurately 
account for such near-discontinuities. This is because the basis 
functions for standard spectral analysis: sinusoids, Legendre or 
Leguerre polynomials, Bessel functions, etc., are global, i.e., 
infinite in length, rather than discrete; whereas fractal analysis 
naturally assumes a discrete Chapeau (top hat) basis function over 
the time resolution in question. And this basis function 
inherently spans the ideal range: from the chosen span of the time 
series to the available limit of resolution. Therefore, functions 
involving frequencies shorter than the available resolution are not 
required to portray large jumps or discontinuities. This suggests 
that fractal analysis may be quite useful in analyzing discrete, 
coherent wave trains and intermittent turbulence structures in 
general . 
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II. THEORY 



A. Current ways to distinguish waves from turbulence 

Since linear waves transport no heat, the cross-spectral method 
assumes for waves that the vertical velocity and temperature 
exhibit a small cospectrum and large quadrature spectrum, i.e., 
temperature lags vertical velocity by 90°, while turbulence shows 
a large cospectrum and small quadrature spectrum. But Finnigan 
(1988) noted that near-surface gravity waves are non-linear, so the 
co-spectrum is usually comparable to the quadrature spectrum. 

The spectral gap method assumes a wave/turbulence gap in the 
boundary layer power spectra (Nai-Ping, 1983). Caughey (1977) ties 
the position of this gap to the Brunt-Vaisala frequency (BVF) , 
which should be the highest frequency of gravity wave that the 
atmosphere will support. This frequency is often, but not always, 
lower than turbulence injection scales. Wave phase velocity also 
adds to the apparent frequency. Thus, spectral gaps in wave number 
may disappear in the frequency domain for time series from a fixed 
sensor (Caughey and Readings , 1975) . 

Phase averaging subtracts the wave and mean components from the 
signal, leaving turbulence as the remainder (Finnigan, ibid). But 
one must assume a wave exists and its frequency must accurately 
estimated prior to the analysis. A surface microbarograph array is 
usually needed to determine wave frequency reliably, but most tower 
sites lack microbarographs. Problems also arise, if wave 
amplitudes change greatly with time, if coherence is lost after a 
few periods, or if non-linear dispersion changes the wave 
frequency. We suggest that phase averaging is suitable for some 
aspects of wave/ turbulence analysis, but only after wave presence 
has been established, perhaps most easily by multi-fractal 
analysis . 

The following discusses how multi-fractal analysis can establish 
such wave presence from single point measurements, without 
microbarographs, or assumptions about quadrature spectra or 
spectral gaps. 

B. Basic concepts of self-similar fractal dimnesion 

For clarity we outline basic fractal concepts before describing the 
self-affine multi-fractal operator used for this analysis. 

The minimum number of boxes of side e needed to cover a set of 
points on a two-dimensional plot scales as 

Nie) = , (II. 1) 

where D is defined as the capacity dimension of the set and F, the 
lacunity (Mandelbrot, 1977) . If the set forms a straight line, 
then D = 1 because, if the box length is halved, twice as many 
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boxes are needed to cover the line. If the set consists of 
uniformly distributed points in a plane, then D = 2, since four 
times as many boxes are needed. 

If the points are not uniformly distributed, then D can be non- 
integer. The set of points is then fractal with dimension D, 
defined in the limit as e approaches zero. This definition is 
called the box dimension, Dg, (Mandelbrot, 1985) , and is given as 



Db 



lim 

t-O 



log N{t) 
logi 



(II. 2) 
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Figure 1 



- Construction of the Cantor set. 
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The Cantor set in fig. 1 has a non-integer Dg. It is formed by 
removing the middle third of a line segment, removing the middle 
third of each remaining segment, and re-iterating to infinity. The 
set is self-similar because the geometric structure seen at the 
largest scale re-appears for all smaller scales removed by 
successive factors of three. For an initial unit length, N(l)=l 
box of side £ = 1 is needed to cover the set. For £=1/3, N(l/3)=2, 
and for £=1/9, N(l/9)=4. Thus, £ is generally 3'", and N(£) is 2", 
where n = 0 , 1 , 2 , 3 , . . . ,oo (Grebogi et al., 1987). Noting that n-*«» 
implies £->0, from eqn. II. 2, Dg for the Cantor set becomes 
intermediate between a point (Dg = 0) and a line (Dg = 1) . That is. 



log S'’ 



^ = 0 .63092 . . . 
log 3 



(II. 3) 



Unlike the analytic Cantor set, Dg must be evaluated numerically 
for most fractal sets. This is done by removing the limit in II. 2 
and reordering terms to get 

logi\/(e) = log F - F^loge . (II.4) 



Then -Dg is the slope of the plot of log N(£) vs. log £, and log F 
is the y-intercept. 

In practice, real time series will contain some noise, so L(£) will 
not be exact. Thus, multiple values of £ should be used, and some 
form of regression employed to find the slope (see section III). 
Moreover, the data may only be fractal over a certain range, so the 
choice of inner and outer scales, £^ and £j, over which the slope is 
determined may be narrower than those defined by purely objective 
limits . 

Digitized data imposes an objective limit on £j because the curve 
shape between discrete data points is unknown, so it is simplest to 
assume a straight line. As before, one plots log N(£) vs log£ to 
find Dg. For £ < £,^ the curve scales as a one dimensional line, 
the assumed shape between adjacent data points. So the inner scale 
of resolution, £j, must span at least three data points, i.e., 2£,„. 
This value happens to be same as the Nykqvist cutoff for Fourier 
spectra, but for different reasons. 

An objective criterion also exists for the choice of outer scale, 
£q, when the time series contains waves. As shown in fig. 2 for a 
single sinusoid, log L(£) versus log £ shows nearly zero slope for 
£ < X/3, whereas for e slightly greater than X/3 the plot is nearly 
vertical. So we suggest for wave containing time series that £<, < 
\^/3, where X^ is the foreseen minimum wave period. 
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Figure 2 - Self-affine L(e) for 1800 point sinusoid. 
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Mandelbrot (1977) describes another way to determine fractal 
dimension for a continuous curve in two dimensions, such as a 
geographical coastline. Here, curve length is defined as 



L(e) = €, N{t) 



(II. 5) 



where N(£) is given as 



Nit) 




(II. 6) 



and F is again the lacunity. 

Lacunity measures the length of one fractally scaling curve to 
another of the same D. I.e., assume that within the limits, to 
£j the fractal dimension is constant for any segment of a real curve 
of finite length. Then for any two segments of unequal length, 
II. 6 shows that the step length ratio will be 



N^ic) _ cD _ (II. 7) 

N^it) 



Using II. 6 in II. 5, L(£) can be written as 

Lit) ~ Ft^-^ . 



(II. 8) 



Taking the log of II. 8 yields the straight line equation. 



logL(e) = logF + (l-D)loge 



(II. 9) 



The slope of log L(£) by log £ yields (1 - D) . To find D first 
find L(£q) . Then, from II. 9, log L(£q) ■= log F + (1-D) log£g. Now 
with L(£j) , log L(£,) = log F + (l-D)log£i. This yields two 
equations in two unknowns, D and F, so in principle 



D = 



log[ 



Lito) 

Lit^) 



] 



log[-^] 



(II. 10) 



This procedure amounts to spanning the curve with circles of 
diameter £. So L(£) is also the minimum number of circles of 
diameter £ required to cover the curve times their diameter. Thus, 
L(£) can be measured by stepping along the curve with a compass 
whose span is set to £. Mandelbrot (1985) terms this estimate the 
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"compass dimension", D^. 

To show that is equivalent to Dg, solve II. 8 for D to get 



D = log I<(e) - log (F) ^ 

log — 
c 

Substituting L(£) from II. 5 gives 

^ logNjt) - log (F) 

logi: 

e 



(II. 11) 



( 11 . 12 ) 



and in the limit as £ 
yields 



0, since N(£) is much larger than F, this 



D = 



lim 

£-0 



log N{t) 

log- 

e 



(11.13) 



which is Dg. 

In essence fractal dimension is a measure of the jaggedness or 
degree of convolution of a curve given in terms of the power law 
dependence of the curve length as a function of resolution. 

The box method divides the data plane into discrete boxes of 
length, £, where N(£) is the number of boxes having at least one 
point on the curve. So, for Dg every box in the domain must be 
checked, even though most boxes are empty, while D(- considers only 
the circles covering points on the curve. So though and Dg are 
identical, Dq is faster to compute. 

C. Time series applications 

Unlike coastlines, a time series plot has axes with different units 
of measure, time versus amplitude. This poses a problem, since the 
"length" of a time series trace cannot be measured without 
specifying an arbitrary scaling ratio between axis units. I.e., 
different ratios will give different values of D^. Thus, though 
calculable, Mandelbrot (1985) suggests that may not be 

meaningful for time series analysis, since its value can even 
exceed two, if the y to x-axis scaling ratio is large enough. This 
warrants exploring a shift in approach to measuring fractal 
dimension, which we describe as follows. 

McHardy and Czerny (1987) applied a "self-affine", "multi-fractal" 
operator to their time series data. Meaning "in the same way", 
self-affinity requires that one axis be stretched or compressed at 
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each successive resolution by the same constant factor with respect 
to the other axis. Thus, McHardy and Czerny defined their self- 
affine "length metric" as 



L(e) = l|^|F(t + e) - F( t) | dt , 
0 



(11.14) 



where now the abscissa scaling changes by the factor, £n/^n+i' each 
time £ changes by the factor £n+i/^n (their length metric differs 
from the standard self-affine length metric, which uses the 
integrand, { [ F ( t+£ ) -F ( t) ] ■ + £-) ''^ ). This new metric leads to, 

^ _ d logL(c) (11.15) 
dloge 



L(£) from 11.14 will be much longer at small £, but note that D in 
11.15 is defined as the change rate of log length with log 
resolution, not the ratio of log length to log resolution. Since 
time, £, cancels in the length metric, L(£) only has units of 
amplitude. This avoids arbitrary ratios between amplitude and time 
units and makes the problem one dimensional, so D is less than 
unity. Thus, eqn. 11.14 presents a natural choice for evaluating 
time series. Hereafter, we refer to D from this method as D^. 

L(£) calculated from 11.14 will contain some noise contribution, 
which is statistically independent of the signal. If its level is 
known, we can estimate its contribution to L(£) by the formula, 
L^bservcd = + L-„„ise • White noise has 0^=1, so it tends to 
increase the fractal dimension of the time series. 

D. Relevant gravity wave and turbulence characteristics. 

Gravity waves in the stable boundary layer can be generated by a 
number of mechanisms, among them: wind shear (Kelvin-Helmholtz 
instability), storm impulses, and flow over an obstacle. Their 
amplitudes can vary from a few centimeters to hundreds of meters, 
with periods from less than a minute up to 40 minutes (Stull, 
1988) . Wave generation due to flow over an obstacle is significant 
here. Hunt (1980) gives the natural wavelength of flow over a hill 
as 



A = 2nUjBVY , (11.16) 

where Uq is mean wind speed. The hill "length" is given as Lj , the 
horizontal distance across the hill at half its maximum elevation. 
If X - 5L] then lee waves may occur with X. Moreover, if £ - 2L, , 
then strong lee waves are possible. Observational studies such as 
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Caughey and Readings (ibid), Caughey (ibid), and Nai-Ping (ibid), 
suggest that waves and turbulence commonly coexist in the nocturnal 
boundary layer, and turbulence generation can sometimes depend on 
the wave. The idea that gravity waves transfer energy to Kelvin- 
Helmholtz waves which then "break" is widely used as a model for 
turbulence generation by waves (Stull, ibid; Atlas et al., 1970), 
though this model is not universally accepted (Hines, 1988) . 

Finnigan (ibid) resolved the wave and turbulent portions of the 
total kinetic energy by phase averaging, and showed horizontal TKE 
developing in the first quarter wave cycle, while vertical TKE 
arose during the ensuing third of the cycle. Both components 
transferred energy from the wave to turbulence. They noted that 
the wave seemed to modulate the turbulence, but the presence of the 
turbulence had no apparent effect on the wave. 

Gossard et al. (1985) suggest a generation mechanism for boundary 
layer turbulence where local gradients of 0, u, and v increase 
steadily, while turbulence decreases. Vertical shear eventually 
reaches an insupportable value of local Richardson number, leading 
to growth of local Kelvin-Helmholtz instabilities which grow into 
turbulence. These various proposals for turbulence initiation and 
wave/turbulence energy transfer are discussed in the context of our 
fractal analysis in section IV. 
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III. METHODS 



A. Data 

The Boulder Atmospheric Observatory (BAO) (Kaimal and Gaynor, 1983) 
supplied the data for this study. Data samples included eight 
levels at; 10, 22, 50, 100, 150, 200, 250, and 300 meters, of: 

1. u, V, and w velocity components sampled at 10 Hz by 
sonic anemometers. 

2. wind direction and speed at 1 Hz by propeller vane 
anemometers . 

3. temperature at 10 Hz by platinum wire thermometers, and 
at 1 Hz by quartz thermometers. 

The 1 Hz data were given only as 10 second averages. The 10 Hz 
data were given at both 10 Hz and as 10 second averages. Noise 
levels were < .01 °C for the platinum wire thermometer, and < .03 
m s'* for the sonic anemometers (J. Gaynor, personal communication) . 

The data and FORTRAN input read programs were sent on 9 track tape 
by John Gaynor and Dave Welch of the NOAA Wave Propagation 
Laboratory. The programs were run on a NPS Computer Science 
Department Sun4 workstation, with analysis on the NPS mainframe. 
The data covered three (MST) time periods: A) 2340 9/7/83 - 0340 
9/8/83; B) 0440 - 0620 9/9/83; and C) 0000 - 0620 9/19/86. 

The temperature time series seemed the least noisy of the available 
data. So 10 Hz temperature data were plotted at levels four and 
five (100 and 150 meters) for all three periods. A "wave" was 
apparent on visual scanning period A at level four, between 0040 
and 0105 MST, September 8, 1983 (see fig. 3, windows 21-27). 
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Figure 3 - Temperature (»C) vs. three minute windows. 
The wave occurs in windows 21 - 27; turbulence episodes 
in windows 45 - 47, 66 - 67, and 73 - 75. 
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since winds were southerly throughout, the "wave" was consistent 
with the above lee wave criteria and flow over a small hill south 
of the tower. I.e., eqn. 11.16 gave X = 1900 meters, for an 
observed BVF = 0.03 Hz, and = 9 m s'’ . From a USGS map for the 
120 ft hill (elev. 5280 ft MSL) , Lj = 300 meters along the tower 
direction; so 5L| ~ 1500 meters. This fits the condition X ~ 5Lj. 
But before and after the "wave" period, the BVF dropped below 0.01 
Hz, so X exceeded 6200 meters, well beyond the lee wave condition. 
This account is consistent with terrain induction of the wave. Lee 
wave conditions were approached for windows 43-53, and 61-66, but 
not as closely, and waves were also not evident in the time series. 
Other periods showed no clear wave evidence. Period C also showed 
pronounced noise spikes at 20 minute intervals and high noise 
levels throughout most of its length. 

In period A a reduced temperature variance region followed the 
"wave" episode. This let us compare fractal characteristics of 
both "wave" and non-wave portions of the time series. Thus, we 
focused on level four data from period A for this study. 

B. Analysis measures 

D^, D^, turbulent kinetic energy (TKE) and its component velocity 
variances, bulk Richardson number (Rg) , Brunt-Vaisala frequency 
(BVF) , buoyancy length (Ig) , and fast fourier transform (FFT) 
spectra were all tested as analysis tools on the u, v, w, and T 
time series at level four. Rg, BVF, and Ig are bulk measures, and 
were taken across levels three to five, a vertical span of 100 
meters. The above parameters were checked for possible correlation 
with D^. 

The main disadvantage of bulk measures is that they apply to the 
entire layer rather than a local sensor; D^, TKE, FFT spectra, 

and variance are taken from local values valid only at the sensor. 
So if an eddy is smaller than the 100 meter layer thickness, it may 
impact the local sensor but not the heights at which the bulk 
gradients of U, V, and 6 are measured. 

was computed for the level 4, 10 Hz temperatures for adjacent 

three minute windows throughout period A. Three minute windows 
were chosen, since the "wave" had a period of about 200 seconds. 
As discussed above, £j was set to three data points (.2 sec). 
was set to the window length, 1800 data points, for the first run. 
For all runs the time axis was expanded by a factor of 100, so that 
.1 seconds on the ordinate was equal to .01°C on the abscissa. 

The first run for showed that the curve scaled linearly (slope 
on the log-log plots was zero) for all log e > 2.25, or e > 17.8 
seconds for all windows. For most windows, the linearly scaling 
portion began above £ > 3.2 seconds. For better resolution on the 
fractal part of the curve, the run was repeated with an = 17.8 

seconds. Figure 4 shows a typical plot from this run. 
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The steep, linear part of the curve was linearly regressed to get 
Dc. The linear region varied from window to window, depending on 
the outer scale. was determined manually for each window, 

according to the following procedure. 

1. A line was drawn with a straight edge along the steep, 
linear portion of the curve (point A in Fig. 5) . 

2. A second line was drawn along the shallow, trailing edge 
of the curve at high values of log e. 

3. The ordinate where these two lines intersected was taken 
to be the outer scale, 

Some windows also required that £j, be determined manually in a 
similar fashion. 
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was also computed for the temperature and the three velocity 
components, using the same three minute windows for Period A. £j 
and £p were again .2 and 180 seconds for the first run. This 
allowed approximately three decades of dynamic range. 

The first run showed for all windows that log L(£) vs. log £ had 

constant, non-zero slope for £ from .2 to around 60 seconds, but 

not beyond this time duration. Thus, only log £ < 2.75 {e ~ 600) 
was used, since this still provided a satisfactory 2.5 decades of 
dynamic range. As mentioned above, another reason for setting 
= 600 is that X„,/3 = 600 for the apparent "wave” period. Thus, 
unlike Dc objective criteria could be used to determine the inner 

and outer scale cutoffs for the slopes. Again a simple linear 

regression was used to get for each window. 

Mean turbulent kinetic energy (TKE) is defined as + <^w") / 

where a^, a^, and are the respective standard deviations of u, 

V, and w. Variances are often measured over 30 minute to one hour 
periods because the apparent "spectral gap" often occurs at about 
1 hour and allows a convenient separation between large and small 
scales. Statistical confidence in the value of second moments such 
as variance also may not reach appropriate levels without averaging 
over tens of minutes. In our case such a long averaging time would 
not be useful in determining small scale differences between waves 
and true turbulence. Instead, the TKE and variances were 
calculated over the same three minute intervals as and for 
purposes of direct comparison. All other parameters such as bulk 
Richardson number or Brunt-Vaisala frequency were also calculated 
over three minute windows. 

Bulk Richardson number, R 3 , was checked for possible correlations 
with D(; or D^. Ideally, Rf, the flux Richardson number, or R,, the 
gradient Richardson number, should be used, since our interest is 
in the local stability at the sensor. Both Rf and Rj require that 
the local vertical gradient of the mean wind be known, but this was 
not available, since the sensors were spaced 50 meters apart 
vertically. So, instead we used the bulk Richardson number, 

g Az 

~ [ (A77)2 + (Av)2] 



where A represents the quantity difference between the bottom and 
top of the layer. The layer thickness, Az , was 100 meters, with 
the layer centered on the level of interest. 

The virtual potential temperature, 6 ^, was assumed equal to 
potential temperature, 0. Though Stull (ibid) emphasizes that 0^ 
can differ from 0 by up to 4°C, this will not seriously affect Rg 
for three reasons. First, d9^/dz will not differ much by 
substituting 0 for 0^; second, a 4°C difference in the denominator 
will not change Rg by more than 10 percent; and third, the absolute 
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value of Rg is not important. In this study Rg is only 
qualitatively important as an indicator of stability changes. 

The Brunt-Vaisala frequency, BVF, was chosen to measure static 
stability and, like Rg, possible correlations with or were 
tested. BVF is defined as 



BVF 




dz 



As with Rg, we assumed that dB^,/dz — AG/Az because we wished to 
explore correlations in temporal changes in BVF with other 
parameters and were not interested in the actual BVF values. 

Ig, the buoyancy length scale, was also determined as the standard 
deviation of the vertical velocity divided by the Brunt-Vaisala 
frequency ( Ig = a^/BVF) , and was assumed to be a measure of the 
dominant eddy scale. 

Detrended FFT spectra were also computed using a cosine squared, or 
Hamming window to excise ringing and high frequency noise. 
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IV. RESULTS AND DISCUSSION 



A. Comparisons with local measures 

Figure 3, the temperature series, shows both "turbulence" and 
"wave" episodes. Using three minute windows, Figure 6a shows 
distinct temperature variance, a-f , peaks. The largest occurs in 
windows 21-28, cotaneous with the "wave". Three other window 
periods: 45-47, 64-67, and 73-75 are cotaneous with "turbulence". 
The other peaks, like window 62, were cotaneous with monotonic 
temperature changes across the window which did not seem to be 
associated with the "wave" or "turbulence". These are absent in 
the and plots (figs. 6b, c) , but do appear in (fig. 6d) . 
Figure 7 shows a TKE spike during the "wave", but also for the 
three "turbulent" episodes. So the TKE and variances do not seem 
to distinguish between "wave" and "turbulence"; all appear as local 
maxima . 

For fig. 8 shows a typical log L(s) versus log £ plot for a 
three minute temperature window. Most plots were quite linear for 
2.5 decades on the ordinate, from £ = .2 to 60 sec. Figure 9 shows 
the temperature Unlike variance, the lowest values occur in 
windows 21-28 during the "wave". Three other local minima occur in 
windows: 42, 63, and 78, cotaneous with periods of seemingly low 
"turbulence" in the time series. Unlike the extended "wave" 
period, the other local minima occurred only in isolated individual 
windows . 

Figure 10 compares the temperature and velocity D^. The curves 
agree, with near-simultaneous minima and maxima on all plots. 
Figure 11 shows scatter plots of from one time series versus 
another, indicating good linear correlation. 
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Figure 7 - Mean turbulent kinetic energy (TKE) . 
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Figure 8 - Typical log-log plot from self-affine 
algorithm (D^^) from the temperature time series. 
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Figure 9 - from temperature time series. 



24 



data) Da (u — data) Da 



o 




0.0 “ I J I I i I I I • ) I I > I ) ) H I 1 1 I i I I ^ ^ I I [ I 1 1 I 1 1 I M I I I I I I I I I I I I I I *rrin' i 't i 1 1 i 1 1 1 1 1 i 1 1 i 1 1 1 1 1 i 1 1 



0 10 20 30 40 50 60 70 80 




Figure 10 - from: (a) T; (b) w; (c) u; and (d) v time 

series. 



25 



1.0 n 




5 ^ 

'^ 0.4 - 



o : 

Q : 

0.2 - 



0.0 



: (a) 



Correlation coefficient = 0.71 



rmTm 1 1 1 rn n n i i rn rrrT r rprr m m i i n m h n i t p m m m n i h rni r rm r n i m m m n t i 



0.1 



0.2 



0.3 



0.4 



0.5 



0.6 



0.7 



Da (T - data) 



0.8 



0.9 



1.0 




0.2 - 





0.0 



(b) 



0.2 



Correlation coefficient = 0.63 

I M 1 1 1 1 1 1 ri 1 1 1 1 1 1 1 1 n n I H I M 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 I T| 

0.3 0.4 0.5 0.6 0.7 0.8 0.9 I.O 

Da (u — data) 



Figure 11 - Scattergrams of from: (a) T vs. w; 

vs. w. 



(b) u 



26 



Thus, seems to discriminate between "turbulence" (high D^) and 
"wave" (low D;^) episodes. For nocturnal boundary layers this 
points to the possible use of ~ 0.35 from temperature, or 
a 0.5 from vertical velocity data as conditional sampling cut- 
off values to distinguish wave dominated periods from time series 
of hot-wire and sonic anemometer data, while retaining most of the 
real turbulence. However, more data containing both waves and 
turbulence must be analyzed before a firm value can be suggested. 
Different cut-off values may be appropriate for time series from 
other instruments, such as hot film, lidar, and cup anemometers 
having different amplitude ranges and time resolutions. 

Though the following suggestion is untried as yet, note that is 
large for "turbulence" but small for "waves", while the variances 
are large for both. So wave periods can be highlighted by scanning 
cotaneous windows with the combined operator, while 

turbulence may be highlighted with the operator, 

All values increase toward the end of the "wave", cotaneous with 
larger high frequency fluctuations in the T and w time series. 
This may indicate "wave-break", or wave-turbulence energy transfer. 
So may also be useful in determining wave-break episodes, unlike 
other purely local measures such as variances or TKE. Waves may 
augment the local shear, thereby instigating instability, but they 
may also enhance TKE by increasing the buoyancy flux. Figure 12 
shows a large increase in the non-phased averaged temperature flux 
and hence buoyancy flux during the latter part of the "wave" 
period, without a concommitant increase in the Reynolds stress 
portion of the TKE shear generation rate. Rg also plummets during 
"wave-break", and before one "turbulence" burst, but not prior to 
the other two "turbulence" episodes. 

Cotaneous with the calm preceding these "turbulent" bursts, two of 
the local minima occur just before the "turbulence" bursts and 
maxima at windows 46 and 66. This is consistent with the 
proposal that local shear may gradually increase during periods of 
relative calm; but when such shear reaches an insupportable level, 
K-H instability may be initiated and rapidly grow into turbulence. 
This shear based initiation differs from the buoyancy based "wave" 
case above. 

The highest peak at window 36 bears more scrutiny, since it is 
cotaneous with a relatively calm appearing, low variance portion of 
the time series following the "wave" period. One possibility is 
that "wave" initiated turbulent mixing may persist for some time as 
intense but small scale turbulence, difficult to resolve by visual 
scanning alone. may be a more acute discriminant than the eye 

itself . 
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Fig. 12 - D^(T), stress, and temperature flux in Period A. 
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Alternatively, one might also suggest, due to the 1/2 weighting, 
that is largely a scale separator which seems to distinguish 
between waves and turbulence only because turbulence in the stable 
boundary layer is generally smaller scale than waves. However, 
this does not explain the minima also seen at windows: 5, 42, 

56, 63, and 78, where the fluctuations seem to be small scale. 

Since Dc, the self-similar, fractal dimension depends on the y to 
x-axis scaling ratio, we arbitrarily set .1 second = 100°C. As 
expected, behaved differently from D^, rising dramatically 

during the "wave" episode, and falling immediately after. This 
same behavior was observed during the "turbulent" episodes; thus, 
D(- did not distinguish between the "wave" and turbulence. Moreover, 
the long, ad hoc algorithm involved a dynamic range, C; - 
smaller and more subjective than the limits used for D^. For these 
reasons does not seem useful as a wave/ turbulence discriminant 
for time series and was not pursued further. 

Figure 13a, b show Fourier spectra of the temperature data before 
and during the wave. Before taking the FFT, the data were linearly 
detrended and tapered. Despite such processing, the plots are 
quite noisy, showing only slight evidence of a spectral gap between 
the "wave" and higher frequency fluctuations. There is a decided 
increase in slope for the lower frequency components above 0.1 Hz 
which may indicate wave presence. However, the "turbulence" periods 
do not exhibit characteristics of any distinctive sort in the FFT 
plots, as shown in fig. 14a, b. 

also has a natural advantage over traditional spectral methods 
in that it assumes that the amplitude is constant over the time 
resolution in question. So actually relies on local Chapeau 
basis functions rather than global basis functions, such as Fourier 
series, Bessel functions, etc. to obtain digitized amplitude 
differences. Thus, as illustrated in fig. 15, has no problem in 
tracking large, nearly discontinuous jumps in amplitude, while FFTs 
require that linear combinations of wave numbers much higher than 
£j be found to both simulate such jumps and suppress fictitious 
ringing outside the area of the jump discontinuity. 
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Figure 13 - FFT spectra from temperature time series: (a) before 
wave; and (b) during wave. 
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Figure 14 - FFT spectra from temperature time series: (a) before 

third turbulent episode; and (b) during turbulence. 
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B. Comparisons with bulk measures 

Figures 16(a) and (b) show the temperature compared with Brunt- 
Vaisala frequency (BVF) . The BVF is highest during the "wave" when 
is small. The three minima at windows 42, 63, and 78 occur 
when the BVF is at or near a local maxima, consistent with the fact 
that a higher BVF (more stable) will tend to dampen turbulence. 
Also, most of the maxima occur when the BVF is either decreasing 
or at a minimum, which indicates a less stable, more turbulent 
atmosphere. But high BVF values occurred during both "wave" and 
"turbulence" periods, so BVF cannot be used to resolve the two. 
Figure 16(c) shows a scatter plot of the temperature versus BVF, 
which indicates a weak negative correlation. 

Figure 17 shows another weak, but positive, correlation between 
and Ig, the buoyancy length scale. The correlation is highest 
during the "wave" episode, when Ig drops rapidly, simultaneous with 
a rapid drop in D^. After the "wave", and Ig both rise. The u 
and T variances and BVF show local maxima at window 62, while the 
u and w variances and Ig show local minima. The temperature time 
series shows a rapid monotonic temperature fall. Ig does not seem 
to distinguish between "wave" and turbulence. All such episodes 
show moderate Ig on the order of 8-12 meters. Since the BVF is in 
the denominator of Ig, Ig peaks during windows 18-20, and 35-40 when 
the BVF is low. This seems to occur when the temperature time 
series is relatively calm. Thus, large Ig does not seem to 
correlate with turbulence, and does not seem to be a good measure 
of dominant eddy scale, as some have suggested (Stull, ibid). 

Comparison of with Rg in fig. 18 shows little correlation except 
during the "wave" when Rg rises dramatically, indicating increased 
stability. None of the bulk measures seem to reliably distinguish 
between "wave" and "turbulence" episodes. Note however, that high 
correlations between and BVF, Rg, or Ig may not be reasonable, 
since is a local measure, whereas the latter three are bulk 

measurements from widely vertically separated sensors. Higher 
correlations might be expected, if more local measurements of Ad 
and Au were available. However, this capability is not currently 
present on the BAG tower or similar facilities. 

The wave/turbulence energy transfer during windows 25-28 may 
initially blend the potential temperatures across levels 3 to 5 
(see Fig. 19), resulting in high Ig and low BVF after window 35. 
This is consistent with the local maxima in the w and v variances 
around window 35. However, again the large maxima in window 35 
of the u, V, w, and T traces may suggest continued turbulent 
blending, principally at smaller scales, even after the level 4 and 
5 potential temperatures have become well-mixed. 

Though not suggested by our data, it is still possible that such 
smaller scale blending is unduly accentuated, since is 

normalized by 1/fi", where n was assumed to equal unity. For 
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inertial eddy cascade processes, n may indeed be less than unity. 
As a first guess, since the turbulence energy spectrum tends to 
have a -2/3 slope in the inertial subrange, one might expect the 
turbulence velocity spectrum to have a -1/3 slope. In that case, 
perhaps the correct self-affinity factor for inertial subrange 
turbulence should be n = 1/3. So to modify one would first have 
to modify the integrand in 11.14 somehow to retain the units of 
pure amplitude for L(e) . 

Secondly, in practice the inertial subrange scaling assumption is 
not entirely viable. The BAO data resolution is 0.1 seconds, which 
enters the dissipation range where the spectral energy slopes are 
much higher. To eliminate contributions from this range, the inner 
scale cutoff would have to be at longer times, thus narrowing the 
dynamic range over which is computed. This may be possible for 
convective data, since a longer window must be used to capture 
fluctuations due to boundary layer scale eddies. 

Thirdly, stability also tends to compress the vertical extent of an 
eddy. Thus, inertial subrange turbulence may not have a -% slope 
in the energy spectrum, as in self-similar turbulence. That is, 
there will be less energy in the larger scales. Mahrt and Gamage 
(1987) have studied vertical/horizontal velocity aspect ratios for 
structure functions as a function of stability. However, their 
stability categories were not well defined by Rj or other measures. 

In light of these difficulties and lack of supporting data, we 
retain the present self-affine operator, Dys^, and note that it 
appears to be the first potentially operationally useful wave/ 
turbulence discriminant. It may also be useful as a fast pre- 
processor to quickly locate intermittent wave or turbulence periods 
for more detailed phase averaged studies. 

As an extension of this study, another potential use for is as 
a general measure of the degree of chaos, applied to systems 
ranging from pure waves and limit cycles, through standard chaotic 
systems such as the Henon, Poincare, and Lorenz maps, as well as 
transitional turbulence, and turbulence displaying extended 
inertial subranges. One effort has begun to test Dy^ on a set of 
differential equations which describe coupled harmonic oscillators. 
By eye the degree of chaos expressed in each oscillator's motion 
seems to increase or decrease monotonically along the oscillator 
lattice. An adequate measure this chaos has not been shown as yet. 
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Figure 16 - (a) BVF; (b) from T time series; (c) 
vs. BVF. 
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IV. SUMMARY AND CONCLUSIONS 



Time series of T, u, v, and w from a stable atmospheric boundary 
layer were examined using fractal techniques. The time series were 
taken from lOHz sonic anemometer and hot wire data at a height of 
100 meters on the mast at the Boulder Atmospheric Observatory 
(BAO) . The time series featured a probably terrain induced gravity 
"wave" period, as well as several periods of intermittent 
turbulence. 

Both self-similar fractal dimension, D^, and self-affine fractal 
dimension, D^, were explored in the analysis. The latter was 
clearly superior in practice as well as principle. was the only 
measure among the eight (five local and three bulk) studied which 
seemed to distinguish the "wave" from "turbulent" episodes, 
reaching minima during "wave" and maxima during "turbulence". 

Since is a local measure, this implies that tower data from a 
single fast response sensor can be used for wave/ turbulence 
discrimination, while phase averaging requires extensive numerical 
processing and microbarograph array data as well. Most tower sites 
are not equipped with microbarograph arrays. 

Phase averaging also requires that the wave period be known prior 
to analysis, otherwise the wave and turbulent components cannot be 
distinguished. also requires an estimate of the minimum 
expected wave period, X,,,, to specify the outer scale for the 
slope. However, the required precision is perhaps a factor of two, 
much lower than is needed for phase averaging. For we note that 
nocturnal atmospheric boundary layer gravity waves have typical 
periods in the range 100 - 300 seconds. 

From the limited results we tentatively suggest a conditional 
sampling cutoff value, ~ 0.35, to remove wave data from a hot 
wire temperature time series. However, more data needs to be 
analyzed to establish a firmer value. Since peaks during both 
waves and turbulence, dividing by should highlight the wave 
periods. Similarly, turbulence should be highlighted by multiplying 
Da by 

Correlations between and various bulk measures of stability were 
investigated, with weak correlations (r < .41) found with Brunt- 
Vaisala frequency (BVF) and buoyancy length (1 b)/ little 
correlation with bulk Richardson number (Rg) . The issue of 
inherently low correlations between local and bulk measures was 
discussed . 

Other possible parameters for distinguishing waves and turbulence 
were investigated. Turbulent kinetic energy (TKE) proved less than 
useful, since it included wave energy. Fast Fourier Transform 
spectra were able to resolve the "wave" event as a change in slope 
at low frequency, but no distinguishing features were found for the 
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"turbulence" episodes. 



In terms of computational efficiency does not require that 

coefficients be determined for linear combinations of wave 
frequencies to approximate moments of the time series amplitude, 
i.e., no transform is needed. Moreover, since the slopes used to 
determine seem well-behaved, perhaps sufficiently accurate 

values of may be determined by computing L(e) only at the inner 
and outer time resolutions. This would render computations 
orders of magnitude faster than FFTs. Moreover, FFTs tend to treat 
finite wave trains, i.e., discontinuities quite poorly. avoids 

this problem because it involves local Chapeau basis functions, 
rather than global transcendental or other basis functions. This 
makes a potentially suitable analysis tool for turbulence 

intermittency and coherent structures in general. 

also shows promise in resolving wave-turbulence energy 
transfers, or "wave-break" events. We suspect that "turbulence" 
following the "wave" period studied may have been initiated by 
growing positive buoyancy flux, while other "turbulence" episodes 
may have relied more on shear generation of TKE. However, 
correlations between and "wave-break" should be tested in more 
detail, perhaps with the Gossard (ibid) data set which involved 
vertical carriage traverses along the BAO tower, wherein local 
values of R; were measured. Ultimately, one would like to establish 
as a function of local stability, so that turbulence cut-off 
values could be associated with critical Rj. 

Tentatively, we interpret the anomalous temperature peak 

following the "wave" period as indicative of intense but small 
scale turbulence, not easily resolved by visually scanning the time 
series trace. Alternatively, due to the 1/e self-affine 
normalization, may act as simply a "scale" discriminator rather 
than a wave/ turbulence discriminator. However, minima were also 
found during periods where visual scanning suggested low turbulence 
levels, not consistent with the latter interpretation. Further 
investigation to test on convective data is suggested to explore 
this possibility. Proceeding in the other direction as well, 
should be tested on systems which lie between the wave and limit 
cycle regimes and turbulent systems displaying extended inertial 
subranges, such as the Henon, Poincare, Lorenz, and other standard 
chaotic attractors, to see if can distinguish degrees of chaos 
in general. An immediate effort is underway to test on 

differential equations which describe coupled harmonic oscillators. 
To the eye in certain modes each oscillator along the lattice seems 
to display increasingly chaotic motion. An adequate measure of the 
this chaos has not been demonstrated as yet. 

An n = 1/3 self-affine scaling factor for £ is discussed as 

possibly more suitable for turbulence, though this would require 
modification of the integrand in 11.14 to avoid making L(£) also a 
function of time. 
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